SARIMA

Ejercicio de Simulación

Vamos a hacer un ejercicio de simulación para ver como se identifica la componente estacional.

library(urca)
library(forecast)
library(tseries)

    ‘tseries’ version: 0.10-51

    ‘tseries’ is a package for time series analysis and
    computational finance.

    See ‘library(help="tseries")’ for details.

library(lmtest)

Loading required package: zoo

Attaching package: ‘zoo’

The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric


Attaching package: ‘lmtest’

The following object is masked _by_ ‘.GlobalEnv’:

    unemployment

library(uroot)
library(fUnitRoots)

Loading required package: timeDate
Loading required package: timeSeries

Attaching package: ‘timeSeries’

The following object is masked from ‘package:zoo’:

    time<-

Loading required package: fBasics

Attaching package: ‘fUnitRoots’

The following objects are masked from ‘package:urca’:

    punitroot, qunitroot, unitrootTable

library(sarima)

Loading required package: stats4

Attaching package: ‘sarima’

The following object is masked from ‘package:stats’:

    spectrum

require("PolynomF")

Loading required package: PolynomF

###Simulación de un proceso con raíz unitaria estacional
#x11()
x <- ts(sarima::sim_sarima(n=144, model = list(iorder=0, siorder=1, nseasons=12, sigma2 = 1),n.start=24),frequency = 12)
plot(x)

acf(x,lag.max = 36)

monthplot(x)

nsdiffs(x)####Decreta cuantas diferencias estacional a través de la aplicación de 

[1] 1

###Algunas pruebas de raíces unitarias estacionales.


###diferencia estacional
Dx=diff(x,lag=12,differences = 1)###lag:periodo s.
plot(Dx)

acf(Dx,lag.max = 36)

monthplot(Dx)

nsdiffs(Dx)

[1] 0

####Simulación de un SAR
#x11()
x1 <- ts(sim_sarima(n=144, model = list(ar=c(rep(0,11),0.8)),n.start=24),frequency=12)
plot(x1)

acf(x1,lag.max = 36)

monthplot(x1)

nsdiffs(x1)

[1] 0

ndiffs(x1)

[1] 0

p <- polynom(c(1,c(rep(0,11),-0.8)))
solve(p)

 [1] -1.0187693+0.0000000i -0.8822801-0.5093846i -0.8822801+0.5093846i
 [4] -0.5093846-0.8822801i -0.5093846+0.8822801i  0.0000000-1.0187693i
 [7]  0.0000000+1.0187693i  0.5093846-0.8822801i  0.5093846+0.8822801i
[10]  0.8822801-0.5093846i  0.8822801+0.5093846i  1.0187693+0.0000000i

abs(solve(p))

 [1] 1.018769 1.018769 1.018769 1.018769 1.018769 1.018769 1.018769
 [8] 1.018769 1.018769 1.018769 1.018769 1.018769

###Note lo cerca que están la raíces de la no estacionariedad del proceso, por eso
####aunque si bien el proceso es estacionario, notamos hay una cercanía a 
####e tener una componente estacional.
####El anterior modelo puede escribirse como:
x2 <- ts(sim_sarima(n=144, model=list(sar=0.8, iorder=0, siorder=0, nseasons=12),n.start=24),frequency = 12)
plot(x2)

acf(x2, lag.max=48)

monthplot(x2)

nsdiffs(x2)

[1] 1

x3=ts(sim_sarima(n=144, model=list(sar=0.7, sma=0.5, iorder=0, siorder=1, nseasons=12),n.start=24),frequency = 12)
plot(x3)

acf(x3, lag.max=48)

monthplot(x3)

nsdiffs(x3)

[1] 1

###diferencia estacional
Dx3=diff(x3,lag=12,differences = 1)###lag:periodo s.
plot(Dx3)

acf(Dx3,lag.max = 36)

monthplot(Dx3)

nsdiffs(Dx3)

[1] 0

#########3
x3=ts(sim_sarima(n=144, model=list(ma=0.6, sma=0.5, iorder=0, siorder=1, nseasons=12),n.start=24),frequency = 12)
plot(x3)

acf(x3, lag.max=48)

monthplot(x3)

nsdiffs(x3)

[1] 1

Dx3=diff(x3,lag=12,differences = 1)###lag:periodo s.
plot(Dx3)

acf(Dx3,lag.max = 36)

monthplot(Dx3)

nsdiffs(Dx3)

[1] 0

Ejemplo Pasajeros

Vamos a ver como se hace el modelamiento completo de la serie de pasajeros.

Iniciaremos con la transformación Box-Cox y las pruebas de raíces Unitarias

[1] 1

Call:
ar(x = lAirpassengers)

Coefficients:
      1        2        3        4        5        6        7        8  
 1.0013  -0.0959   0.0329  -0.0252   0.0231   0.0005  -0.0157  -0.0762  
      9       10       11       12       13  
 0.1071  -0.0236   0.0675   0.4536  -0.4854  

Order selected 13  sigma^2 estimated as  0.01336

    Augmented Dickey-Fuller Test

data:  lAirpassengers
Dickey-Fuller = -2.147, Lag order = 13, p-value = 0.5152
alternative hypothesis: stationary

Warning: p-value greater than printed p-value

Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 12
  STATISTIC:
    Dickey-Fuller: 3.7872
  P VALUE:
    0.99 

Description:
 Thu Oct 27 08:19:54 2022 by user: 


############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression none 


Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.103864 -0.023652 -0.001455  0.022160  0.126649 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
z.lag.1       0.005439   0.001436   3.787 0.000241 ***
z.diff.lag1  -0.198813   0.072080  -2.758 0.006737 ** 
z.diff.lag2  -0.273730   0.073309  -3.734 0.000292 ***
z.diff.lag3  -0.233607   0.072258  -3.233 0.001589 ** 
z.diff.lag4  -0.293133   0.073926  -3.965 0.000126 ***
z.diff.lag5  -0.206562   0.072058  -2.867 0.004915 ** 
z.diff.lag6  -0.266919   0.071493  -3.734 0.000292 ***
z.diff.lag7  -0.234526   0.071847  -3.264 0.001436 ** 
z.diff.lag8  -0.327393   0.073197  -4.473 1.79e-05 ***
z.diff.lag9  -0.198455   0.073623  -2.696 0.008054 ** 
z.diff.lag10 -0.279931   0.072710  -3.850 0.000192 ***
z.diff.lag11 -0.176122   0.073011  -2.412 0.017394 *  
z.diff.lag12  0.627402   0.072683   8.632 3.34e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04308 on 118 degrees of freedom
Multiple R-squared:  0.8564,    Adjusted R-squared:  0.8406 
F-statistic: 54.14 on 13 and 118 DF,  p-value: < 2.2e-16


Value of test-statistic is: 3.7872 

Critical values for test statistics: 
      1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62


Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 12
  STATISTIC:
    Dickey-Fuller: -1.5325
  P VALUE:
    0.7711 

Description:
 Thu Oct 27 08:19:54 2022 by user: 

Call:
ar(x = dlAirpassengers)

Coefficients:
      1        2        3        4        5        6        7        8  
-0.0993  -0.1226  -0.2824  -0.2875  -0.1184  -0.2406  -0.1829  -0.2971  
      9       10       11       12       13       14       15  
-0.1040  -0.2627  -0.1259   0.5734   0.0179  -0.1667   0.1200  

Order selected 15  sigma^2 estimated as  0.002988

    Augmented Dickey-Fuller Test

data:  dlAirpassengers
Dickey-Fuller = -3.3376, Lag order = 15, p-value = 0.0678
alternative hypothesis: stationary


Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 12
  STATISTIC:
    Dickey-Fuller: -1.2228
  P VALUE:
    0.2243 

Description:
 Thu Oct 27 08:20:04 2022 by user: 

Identificación de la componente ARMA estacional y la componente ARMA ordinaria

####################################
library(uroot)
require(forecast)
######Diferencia Estacional(continuación AirPassengers)#######
monthplot(dlAirpassengers)
nsdiffs(dlAirpassengers)

[1] 1

nsdiffs(AirPassengers)

[1] 1

DdlAirpassengers=diff(dlAirpassengers,lag=12)###lag=s
#x11()
par(mfrow=c(2,1))

plot(dlAirpassengers)
plot(DdlAirpassengers)

monthplot(DdlAirpassengers)
nsdiffs(DdlAirpassengers)

[1] 0

##Autocorrelogramas
#x11()
acf(DdlAirpassengers)

acf(DdlAirpassengers,lag.max = 48, ci.type='ma')# q=0,1, Q=0,1
pacf(DdlAirpassengers,lag.max = 48) # p=0,1,2,...,9, P=0,1

#SARIMA(p=0,d=1,q=1)x(P=0,D=1,Q=1)s=12

Ajuste del Modelo y Análisis de Residuales

##Ajuste del modelo
###Arima Estacional o SARIMA(p=0,d=1,q=3)x(P=0,D=1,Q=1)s=12 con transformación logaritmica

#Modelo MA(1) estacional
modelo = Arima(AirPassengers, c(0, 1, 1),seasonal = list(order = c(0, 1, 1), period = 12),lambda = 0)
coeftest(modelo)

z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ma1  -0.401828   0.089644 -4.4825 7.378e-06 ***
sma1 -0.556945   0.073100 -7.6190 2.557e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

modeloalter= Arima(AirPassengers, c(1, 1, 0),seasonal = list(order = c(1, 1, 0), period = 12),lambda = 0)

## Análisis de residuales
#x11()
residuales <- modelo$residuals
plot(residuales)

acf(residuales,lag.max = 24)

pacf(residuales,lag.max = 24)
#Test de autocorrelaci?n
Box.test(residuales, lag = (length(residuales)/4), type = "Ljung-Box", fitdf = 2)

    Box-Ljung test

data:  residuales
X-squared = 37.874, df = 34, p-value = 0.2969

######Análisis de Outliers
#Test de normalidad
jarque.bera.test(residuales)

    Jarque Bera Test

data:  residuales
X-squared = 5.2265, df = 2, p-value = 0.0733

###Estad?ticas CUSUM
res=residuales
cum=cumsum(res)/sd(res)
N=length(res)
cumq=cumsum(res^2)/sum(res^2)
Af=0.948 ###Cuantil del 95% para la estad?stica cusum
co=0.14422####Valor del cuantil aproximado para cusumsq para n/2
LS=Af*sqrt(N)+2*Af*c(1:length(res))/sqrt(N)
LI=-LS
LQS=co+(1:length(res))/N
LQI=-co+(1:length(res))/N
par(mfrow=c(2,1))

plot(cum,type="l",ylim=c(min(LI),max(LS)),xlab="t",ylab="",main="CUSUM")
lines(LS,type="S",col="red")

lines(LI,type="S",col="red")
#CUSUM Square
plot(cumq,type="l",xlab="t",ylab="",main="CUSUMSQ")                      

lines(LQS,type="S",col="red")                                                                           
lines(LQI,type="S",col="red")

Pronóstico

#x11()
Pronosticos=forecast(modelo,h=12,level=0.95)
plot(Pronosticos)
predic<-predict(modelo,n.ahead=12)
plot(predic$pred)


#####Comparación de pronósticos####
library(fpp)
train <- window(AirPassengers,start=c(1949,01),end=c(1959,12))
test <- window(AirPassengers,start=c(1960,01),end=c(1960,12))
fitmodelo <- Arima(AirPassengers, c(0, 1, 1),seasonal = list(order = c(0, 1, 1), period = 12),lambda = 0)
refit <- Arima(AirPassengers, model=fitmodelo)
fc <- window(fitted(refit), start=c(1960,1))


h <- 1
train <- window(AirPassengers,start=c(1949,01),end=c(1959,12))
test <- window(AirPassengers,start=c(1960,01),end=c(1960,12))
n <- length(test) - h + 1
fitmodelo <- Arima(AirPassengers, c(0, 1, 1),seasonal = list(order = c(0, 1, 1), period = 12),lambda = 0)
fc <- ts(numeric(n), start=c(1960,01), freq=12)
for(i in 1:n)
{  
  x <- window(AirPassengers, end=c(1959, 12+(i-1)))
  refit <- Arima(x, model=fitmodelo)
  fc[i] <- forecast(refit, h=h)$mean[h]
}
dife=(test-fc)^2
ecm=(1/(length(test)))*sum(dife)
ecm

library(tidyverse)
library(tidyquant)
library(TSstudio)
library(SLBDD)
library(timetk)
library(fable)
library(feasts)
AirPassengers_tsb=as_tsibble(AirPassengers)
Ajuste_tsibble_Pass<- AirPassengers_tsb %>%
  model(
    sarima011011_Pass = ARIMA(log(value) ~ pdq(0,1,1) + PDQ(0,1,1)),
    auto_pass = ARIMA(log(value), stepwise = TRUE, approx = FALSE), 
  dummy_Pass=ARIMA(log(value) ~ pdq(0,1,1)+ season()),
 fourier_Pass=ARIMA(log(value) ~ pdq(0,1,1)+ fourier(K = 2))
 )

glance(Ajuste_tsibble_Pass)

#Ajuste_tsibble_Pass %>% select(sarima011011_Pass)%>%
 # gg_tsresiduals()

#augment(Ajuste_tsibble_Pass) %>%
#  filter(.model=='sarima011011_Pass') %>%
#  features(.innov, ljung_box, lag = 10, dof = 3)

Ajuste_tsibble_Pass %>%
  fabletools::forecast(h=12) %>%
  autoplot(AirPassengers_tsb)

Ajuste_tsibble_Pass %>%
  fabletools::forecast(h=12) %>%
  filter(.model=='sarima011011_Pass') %>%
  autoplot(AirPassengers_tsb)

Ajuste_tsibble_Pass %>%
  fabletools::forecast(h=12) %>%
  filter(.model=='dummy_Pass') %>%
  autoplot(AirPassengers_tsb)

Ajuste_tsibble_Pass %>%
  fabletools::forecast(h=12) %>%
  filter(.model=='fourier_Pass') %>%
  autoplot(AirPassengers_tsb)

Vale la pena decir que si los residuales son no correlacionados pero ellos no tiene distribución normal, se puede usar los intervalos bootatrap ver sección 5.5 libro fpp3. Esto se hace añadiendo bootstrap=TRUE en la función forecast().

Intervenciones

Accidentes laborales desde enero de 1979 hasta diciembre de 1998. Aquí vamos a tener en cuenta tres variables de intervención: * \(D_t:\) número de días laborales. * \(I_{1t}^{(h)}:\) variable de impulso que toma el valor 1 si la semana santa cae en ese mes, y cero en otro caso. * \(I_{2t}^{(h)}:\) Variable que toma el valor de 1 cuando febrero tiene un día más por ser un año bisiesto.

La razón por la cual es razonable incluir estas variables es porque el número de accidentes laborales puede estar siendo afectado por estos sucesos. Sin embargo, hay mas formas de incluir tales efectos, por ejemplo la función calendar.effects del paquete tsoutliers extrae estos efectos pero desde otra respectiva. Por ejemplo, para los días laborales se calcula una variable asociada llamada día laboral(\(td_t\)), la cual es como sigue: Sea \(wd_t\) el número de días laborales por ejemplo en un mes t. \(nwd_t\) es el número de días no laborales en el mes \(t\). La variable \(td_t\) es como sigue: \[wd_t=\# \text{días laborales}-festivos\]

\[nwd_t=\# \text{días no laborales}+festivos\] y así

\[td_t=wd_t-(5/2)nwd_t.\] Los días laborales son los días de lunes a viernes, y los días no laborales son los sábados y domingos. El número de días festivos en el tiempo \(t\) se pasan como un objeto para la función a través del argumento holidays. El efecto de semana santa se define como la porción de días antes de la semana santa(de un total de 6 por defecto) que caen en marzo y abril.Los restantes meses se les asigna el valor de cero.

El año bisiesto se asocia con un valor de 0.75 al mes de febrero si el año es bisiesto, en caso contrario es -0.25. Los demás meses se definen con cero.

Accidentes_tsibble%>%autoplot()

Plot variable not specified, automatically selected `.vars = Accidentes`

Acá construimos las variables para los efectos de calendario. Adicionalmente, vamos a incluir una variable cambio de nivel para febrero de 1988, ya que se introdujo una nueva definición de accidente laboral en la ley Española. Esta variable de paso se define como sigue: \[S_t^{(110)}=\begin{cases} 0 & \text{si } t<110,\\ 1 & \text{si } t\geq 110. \end{cases} \] El valor de 110 es porque corresponde a la observación número 110 en el tiempo. Esta variable se puede crear usando la función filter.

Una función cambio transitorio, tambien puede incluirse \[C(t)=\frac{\omega B}{1-\delta B}S_{t}^{(h)}\]

como sigue “tc_0_4 <- filter(tc, filter = 0.4, method =”recursive”)”

Vamos ajustar un modelo basado en el logaritmo del número de accidentes. Un modelo incial propuesto sin las regresoras es: \[SARIMA((2,1,0)\times(1,1,1))\]

forecast::forecast(modelo_accidentes_xreg,h=1,xreg=regresoras_futuras)

Warning: xreg contains different column names from the xreg used in training. Please check that the regressors are in the same order.

ajuste_final<-aus_cafe%>%model(
  `FourierK=1`=ARIMA(log(Turnover)~fourier(K=1)+pdq(1,1,0)+PDQ(0,0,0),),
  `FourierK=2`=ARIMA(log(Turnover)~fourier(K=2)+pdq(1,1,0)+PDQ(0,0,0)),
  `FourierK=3`=ARIMA(log(Turnover)~fourier(K=3)+pdq(1,1,0)+PDQ(0,0,0)),
  `Dummy`=ARIMA(log(Turnover)~season()+pdq(1,1,0)+PDQ(0,0,0))                              )